#ifndef _FUNC_H_
#define _FUNC_H_

#include "Core/Vec.h"
#include <cmath>

#define PI 3.141592653589793238

template <int Dim> class Func;

template <int Dim>
class Func{
public:
  using rVec = Vec<Real,2>;

  virtual const Real operator()(const Vec<Real,Dim>& pt) const = 0;

  virtual const Real F(const Vec<Real,Dim>& pt) const{return 0.0;};

  virtual const Real Fx(const Vec<Real,Dim>& pt) const{return 0.0;};

  virtual const Real Fy(const Vec<Real,Dim>& pt) const{return 0.0;};

  virtual const Real rev1(const Vec<Real,Dim>& pt) const{return 0.0;};

  virtual const Real rev2(const Vec<Real,Dim>& pt) const{return 0.0;};
  
};

class D2Func1: public Func<2>{
public:
  const Real operator()(const rVec& pt) const {
    Real x = pt[0];
    return (cos(x)*cos(x)-sin(x)+1)*F(pt);
  }

  const Real F(const rVec& pt) const {
    Real x = pt[0];
    Real y = pt[1];
    return exp(y+sin(x));
  }

  const Real Fx(const rVec& pt) const {
    Real x = pt[0];
    return cos(x)*F(pt);
  }

  const Real Fy(const rVec& pt) const {
    return F(pt);
  }
};

class D2Func2: public Func<2>{
public:
  const Real operator()(const rVec& pt) const {
    Real x = pt[0];
    Real y = pt[1];
    return (x*x+y*y)*exp(x*y);
  }

  const Real F(const rVec& pt) const {
    Real x = pt[0];
    Real y = pt[1];
    return exp(x*y);
  }

  const Real Fx(const rVec& pt) const {
    Real x = pt[0];
    Real y = pt[1];
    return y*exp(x*y);
  }

  const Real Fy(const rVec& pt) const {
    Real x = pt[0];
    Real y = pt[1];
    return x*exp(x*y);
  }
};

class D2Func3: public Func<2>{
public:
  const Real operator()(const rVec& pt) const {
    Real x = pt[0];
    Real y = pt[1];
    return x+y;
  }

  const Real F(const rVec& pt) const {
    Real x = pt[0];
    Real y = pt[1];
    return x*x*x/6+y*y*y/6;
  }

  const Real Fx(const rVec& pt) const {
    Real x = pt[0];
    return x*x/2;
  }

  const Real Fy(const rVec& pt) const {
    Real y = pt[1];
    return y*y/2;
  }
};

class D1Func1: public Func<1>{
public:
  const Real operator()(const Vec<Real,1>& pt) const{
    Real x = pt[0];
    return -PI*PI*sin(PI*x);
  }

  const Real F(const Vec<Real,1>& pt) const{
    Real x = pt[0];
    return sin(PI*x);
  }

  const Real Fx(const Vec<Real,1>& pt) const{
    Real x = pt[0];
    return PI*cos(PI*x);
  }

  const Real Fy(const Vec<Real,1>& pt) const{
    return 0.0;
  }
  
};

class D1Func2: public Func<1>{
public:
  const Real operator()(const Vec<Real,1>& pt) const{
    Real x = pt[0];
    return exp(x);
  }

  const Real F(const Vec<Real,1>& pt) const{
    Real x = pt[0];
    return exp(x);
  }

  const Real Fx(const Vec<Real,1>& pt) const{
    Real x = pt[0];
    return exp(x);
  }

  const Real Fy(const Vec<Real,1>& pt) const{
    return 0.0;
  }
  
};

class D2Func0: public Func<2>{
public:
  const Real operator()(const rVec& pt) const {
    return 0.0;
  }

  const Real F(const rVec& pt) const {
    return 0.0;
  }

  const Real Fx(const rVec& pt) const {
    return 0.0;
  }

  const Real Fy(const rVec& pt) const {
    return 0.0;
  }
};

class D1Func0: public Func<1>{
public:
  const Real operator()(const Vec<Real,1>& pt) const {
    return 0.0;
  }

  const Real F(const Vec<Real,1>& pt) const {
    return 0.0;
  }

  const Real Fx(const Vec<Real,1>& pt) const {
    return 0.0;
  }

  const Real Fy(const Vec<Real,1>& pt) const {
    return 0.0;
  }
};

class BCFunc: public Func<1>{
public:
  const Real operator()(const Vec<Real,1>& pt) const {
    //return 0.0;
    return (1.0/16.0)*sin(PI*pt[0]);
  }

  const Real F(const Vec<Real,1>& pt) const {
    return 0.0;
  }

  const Real Fx(const Vec<Real,1>& pt) const {
    return (1.0/16.0)*PI*cos(PI*pt[0]);
  }

  const Real Fy(const Vec<Real,1>& pt) const {
    return 0.0;
  }

  const Real rev1(const Vec<Real,1>& pt) const{
    return asin(16.0*pt[0])/PI;
  }

  const Real rev2(const Vec<Real,1>& pt) const{
    return 1-asin(16.0*pt[0])/PI;
  }
};


#endif // _FUNC_H_
